/*
 * MyThread.java
 *
 * Created on November 16, 2004, 2:32 PM
 */

package edu.umd.essic.maskbuilder;


import java.io.*;
import java.awt.geom.Point2D;
import java.awt.Point;
import java.lang.Math;
import java.nio.*;
import java.nio.channels.*;
import javax.swing.JOptionPane;

public class IThread extends Thread{
    Int_Dialog IDialog;
    FileContentsDlg FileDlg;
    private boolean theEnd  = false;
    private boolean Interpolate = true;
    private boolean Splines = true;
    private boolean Bilinear = false;
    private boolean Averaging = false;
    private boolean AdjustMsk = false;
    private boolean Old_Format = false;
    private boolean DoublePrecision = false;
    private static int FUNC, OldFUNC;
    private static String Magic_ID = "DYM2"; 
    private Point2D.Double[][] coord, metric;
    private int[][] msk, mskt;
    private int size_x, size_y, nlon, nlat, nlevel, nlevel_original, nseriestoskip, misval, MVmethod;
    private float startdate, enddate;
    private double minValue, maxValue;
    private float[] xlon, ylat, zlevel, xlon_im, ylat_im;
    private File inputbin, outputbin;
    private int[] xind, yind; // indices of the new grid such that coord[ind]<x[ind]<coord[ind+1]
    private float[][] xcrd, ycrd, var;
    private double[][] vmid, vfin; 
    private float[][][] all_series;
    private double[] y2,ytmp,y2tmp,yytmp,u;  // auxiliary arrays for interpolation
    private int nx_pr, ny_pr;               // for interpolation: size of matrices to be interpolated
    private boolean error = false;
    private boolean paused = false;
    public boolean mask_changed = false;
    public int[][] mask;
    
    /** Creates a new instance of MyThread */
    public IThread(FileContentsDlg FileDlg, File inputBin, File outputBin, int func, int[][] msk, int nx, int ny, Point2D.Double[][] crd, double[][] mx, double[][] my, boolean mask_type) {
       FUNC = func;
       size_x = nx;
       size_y = ny;
       inputbin = inputBin;
       outputbin= outputBin;
       mask = new int[size_x][size_y];
       coord = new Point2D.Double[size_x][size_y];
       metric = new Point2D.Double[size_x][size_y];
       for(int ix=0;ix<size_x;ix++){
            for(int jy=0;jy<size_y;jy++){
               coord[ix][jy]  = new Point2D.Double(0,0);
               metric[ix][jy] = new Point2D.Double(1,1);
           }
        }    
       for(int i=1;i<size_x;i++){
            for(int j=1;j<size_y;j++){
               metric[i][j].x = mx[i][j];
               metric[i][j].y = my[i][j];
           }
       } 
       
       for (int i=0; i<size_x; i++)
           for(int j=0; j<size_y; j++){
               mask[i][j] = msk[i][j];
               //-------------------------------------------
               //Dym files store the coordinates, where data
               //are available, so, we choose centers or nodes
               //according to designed mask_type, true or 
               //false correspondingly
               //-------------------------------------------
               if (mask_type){
                   coord[i][j].x = (crd[i][j].x+crd[i+1][j].x)/2;
                   coord[i][j].y = (crd[i][j].y+crd[i][j+1].y)/2;
               } else {
                   coord[i][j].x = crd[i][j].x;
                   coord[i][j].y = crd[i][j].y;
               }
       }
       if  (FirstFourBytesInt())
           Old_Format = true;
       else
           Old_Format = false;
       if (Old_Format) 
           ReadOldDymFileHeader();
       else
           ReadDymFileHeader();
       var  = new float[nlon][nlat];       
       
       //FileDlg = new FileContentsDlg(null,true);
       //FileDlg.setLocation(20,20);
       IDialog = new Int_Dialog(null);
       IDialog.setLocation(20,20);
       IDialog.jProgressBar1.setStringPainted(true);
       // Initializing dialog contents---------
       FileDlg.nlon_new  = size_x;
       FileDlg.nlat_new  = size_y;
       FileDlg.new_dym   = !Old_Format;
       FileDlg.funcID    = OldFUNC;
       FileDlg.nlon      = nlon;
       FileDlg.nlat      = nlat;
       FileDlg.nlevel    = nlevel;
       FileDlg.minval    = (float)minValue;
       FileDlg.maxval    = (float)maxValue;
       FileDlg.startdate = startdate;
       FileDlg.enddate   = enddate;
       FileDlg.xlon_0    = xlon[0];
       FileDlg.xlon_nlon = xlon[nlon-1];
       FileDlg.ylat_0    = ylat[0];
       FileDlg.ylat_nlat = ylat[nlat-1];
       FileDlg.zlevel_0  = zlevel[0];
       FileDlg.zlevel_nlevel = zlevel[nlevel-1];
       //---------------------------------------
       FileDlg.show();
       if (FileDlg.OnOK){
           FileDlg.getReturnStatus();
           int nlon_return = FileDlg.nlon_new;
           int nlat_return = FileDlg.nlat_new;
           nseriestoskip   = ((int)FileDlg.startdate-(int)startdate)*12;
           startdate = FileDlg.startdate;
           enddate   = FileDlg.enddate;
           zlevel    = null;
           nlevel_original = nlevel;
           if (enddate <= 12){ // i.e. unites are months, for monthly climatology files
               nlevel = (int)enddate-(int)startdate+1;
               zlevel    = new float[nlevel];
               for(int n=0; n<nlevel; n++)
                   zlevel[n] = (float) (startdate+n);
           }
           else {
               nlevel = Math.round(((enddate-startdate)*12)+1);
               zlevel = new float[nlevel];
               for(int n=0; n<nlevel; n++)
                   zlevel[n] = (float) (startdate+n/12.0);
           }
           if ((nlon_return < size_x)||(nlat_return < size_y)){
               FUNC = 0;
               size_x = nlon_return;
               size_y = nlat_return;
               float deltax = (xlon[nlon-1]-xlon[0])/(size_x-1);
               float deltay = (ylat[nlat-1]-ylat[0])/(size_y-1);
               for(int ix=0;ix<size_x;ix++){
                   for(int jy=0;jy<size_y;jy++){
                        coord[ix][jy] = new Point2D.Double(xlon[0]+ix*deltax,ylat[0]+jy*deltay);
                   }
               }
           } 
           if ((nlon_return > size_x)||(nlat_return > size_y)){
               JOptionPane.showMessageDialog(null,"New nlon and nlat should be equal or less then mask size","Dimension error",JOptionPane.ERROR_MESSAGE);
               FileDlg.doClose();
               FileDlg = null;
               error = true;
               return;
           }
           Old_Format = !FileDlg.new_dym;
           Interpolate= FileDlg.interpol;
           Splines = FileDlg.cubicspl;
           Bilinear = FileDlg.bilinear;
           Averaging = FileDlg.averaging;
           AdjustMsk = FileDlg.adjustmsk;
           DoublePrecision = FileDlg.doubles;
           mask_changed = AdjustMsk;
           misval = FileDlg.missingID;
           MVmethod = FileDlg.MVM_flag;
           //ReadGridData("C://Inna//R_work//nav_lon1.txt","C://Inna//R_work//nav_lat1.txt");
           InitInterpolation();
           String title = "Data interpolation: " + inputbin.getName();
           IDialog.setTitle(title);
           IDialog.show();
       } else {theEnd = true;}
       //all_series = new float[size_x][size_y][nbt_total];
       FileDlg.doClose();
       FileDlg = null;
    }
       
    public boolean isError(){
        return error;
    }
    
    public void run() {  
        int t = 0;
        boolean FileMode = true;
        while((IDialog.getReturnStatus()) && (!theEnd)){  
            RunInterpolation(t,FileMode);
            t++;
            FileMode = false;
            int prDone = t*100/nlevel;
            IDialog.jProgressBar1.setValue(prDone);
            if (t == nlevel) theEnd = true;
         }
        IDialog.doClose();
        IDialog = null;
    }
 
    private void RunInterpolation(int t, boolean FileMode){
        float min, max;
        int nbytetoskip = 0;
        if (Old_Format) 
            nbytetoskip = (int)(5+nlat+nlon+nlevel_original+nlat*nlon*(t+1+nseriestoskip))*4;
        else
            nbytetoskip = (int)(9+nlevel_original+(nlat*nlon)*(t+3+nseriestoskip))*4;
            //nbytetoskip = (int)(9+nlevel_original+((nlat+2)*(nlon+2))*(t+3+nseriestoskip))*4;
        ReadDymFile(inputbin,nbytetoskip);
        if (Interpolate){
            if (!Old_Format){// because in the new format mask is not tranposed
                min = Utils.getMIN(var, nlon, nlat, mskt);
                max = Utils.getMAX(var, nlon, nlat, mskt);
            } else {
                min = Utils.getMIN(var, nlon, nlat, msk);
                max = Utils.getMAX(var, nlon, nlat, msk);
            }
            Interpolate(var);
            ScaleData(vfin, size_x, size_y, min, max);
            WriteDymFile(vfin,FileMode);
        } else {
            FUNC = OldFUNC;
            WriteDymFile(var,FileMode);
        }
    }

    private boolean FirstFourBytesInt() {
        int val = 0;
        try{
            FileInputStream fileInput = new FileInputStream(inputbin);
            DataInputStream in = new DataInputStream(fileInput);
            val = Utils.readInt(in);
        }catch (IOException e) {
            Utils.FileReadError(inputbin.getName());
        }
        if ((val<10000)&&(val>0))
            return true;
        else return false;
    }
    
   
    private void ReadDymFileHeader(){
      ByteBuffer dim_buffer = null;
      int buf1len = 9*4; // first read 9 variables 
      dim_buffer = ByteBuffer.allocate(buf1len);
      try {
        FileInputStream fileInput = new FileInputStream(inputbin);
        FileChannel fc = fileInput.getChannel();
        fc.read(dim_buffer); 
        dim_buffer.order(ByteOrder.LITTLE_ENDIAN);
        dim_buffer.rewind();
        
        // === Magic_ID =====================
        int val = dim_buffer.getInt();
        // Magic_ID = String.valueOf(val);

        // === Stretching function ID =======
        OldFUNC = dim_buffer.getInt();

        // === MIN_VALUE ====================
        minValue = dim_buffer.getFloat();
        //maxValue = dim_buffer.getFloat();
        
        // === MAX_VALUE ====================
        maxValue = dim_buffer.getFloat();
        //minValue = dim_buffer.getFloat();
        //minValue = 0.0;
        
        // === NLON =========================
        nlon = dim_buffer.getInt();
        
        // === NLAT =========================
        nlat = dim_buffer.getInt();
        
        // === NLEVEL =======================
        nlevel = dim_buffer.getInt();
        
        // === STARTING DATE ================
        startdate = dim_buffer.getFloat();
        
          // === ENDING DATE ==================
        enddate = dim_buffer.getFloat();
        
        // now we know the size of arrays and
        // the buffer size that we should allocate
        // for the rest of the file header
        xcrd = new float[nlat][nlon];
        ycrd = new float[nlat][nlon];
        msk = new int[nlat][nlon];
        mskt = new int[nlon][nlat];
        xlon = new float[nlon];
        ylat = new float[nlat];
        zlevel = new float[nlevel];
        dim_buffer = null;
        int buf2len = (int) 4 * (3 * nlon * nlat + nlevel);
        dim_buffer = ByteBuffer.allocate(buf2len);
        fc.read(dim_buffer, buf1len);
        dim_buffer.order(ByteOrder.LITTLE_ENDIAN);
        dim_buffer.rewind();
        
        // === X GRID COORDINATES ===========
        for (int i = 0; i < nlat; i++) {
            for (int j = 0; j < nlon; j++) {
              xcrd[i][j] = dim_buffer.getFloat();
              if (xcrd[i][j]<0) 
                  xcrd[i][j] += 360.0;
              if (i == 0)
                xlon[j] = xcrd[0][j];
            }
          }
          
        // === Y GRID COORDINATES ===========
        for (int i = 0; i < nlat; i++) {
            for (int j = 0; j < nlon; j++) {
               ycrd[i][j] = dim_buffer.getFloat();
               //ycrd[i][j] = -dim_buffer.getFloat();
            }
            ylat[i] = ycrd[i][0];
          }
        //=== ZLEVEL =======================
         for (int k=0;k<nlevel;k++){
           zlevel[k] = dim_buffer.getFloat();
         }
         enddate = zlevel[nlevel-1];
        
        //=== MASK =========================
        for (int i=0;i<nlat;i++){
            for (int j=0;j<nlon;j++){
                msk[i][j]= dim_buffer.getInt();
                mskt[j][i]= msk[i][j];
            }
        }
         dim_buffer = null;
      } catch (IOException e) {
          Utils.FileReadError(inputbin.getName());           }
   }
    
   private void ReadOldDymFileHeader(){
      ByteBuffer dim_buffer = null;
      int buf1len = 5*4; // first read 5 variables 
      dim_buffer = ByteBuffer.allocate(buf1len);
      try {
        FileInputStream fileInput = new FileInputStream(inputbin);
        FileChannel fc = fileInput.getChannel();
        fc.read(dim_buffer); 
        dim_buffer.order(ByteOrder.LITTLE_ENDIAN);
        dim_buffer.rewind();
        
         //=== NLON =========================
         nlon = dim_buffer.getInt();
        
        //=== NLAT ========================= 
         nlat = dim_buffer.getInt();
        
        //=== NLEVEL =======================
         nlevel = dim_buffer.getInt();
                  
        //=== STARTING DATE ================
         startdate = dim_buffer.getFloat();
        
        //=== ENDING DATE ==================
         enddate   = dim_buffer.getFloat();
        
         // now let's allocate the memory
         // for file data arrays and
         // for the rest of the file header
         msk  = new int[nlon][nlat];
         xlon = new float[nlon];
         ylat = new float[nlat];
         zlevel = new float[nlevel];
         dim_buffer = null;
         int buf2len = (int) 4*(nlon+nlat+nlon*nlat+nlevel);
         dim_buffer = ByteBuffer.allocate(buf2len);
         fc.read(dim_buffer,buf1len);
         dim_buffer.order(ByteOrder.LITTLE_ENDIAN);
         dim_buffer.rewind();

         //=== X GRID COORDINATES ===========
         for (int i=0;i<nlon;i++){
                xlon[i] = dim_buffer.getFloat();
         }
         
        //=== Y GRID COORDINATES ===========
         for (int j=0;j<nlat;j++){
                ylat[j] = dim_buffer.getFloat();
         }
         //=== ZLEVEL =======================
         for (int k=0;k<nlevel;k++){
              zlevel[k] = dim_buffer.getFloat();
         }
         enddate = zlevel[nlevel-1];
        
        //=== MASK =========================
        for (int j=0;j<nlat;j++){
                for (int i=0;i<nlon;i++){
                    msk[i][j]= dim_buffer.getInt();
                }
        }
         dim_buffer = null;
      } catch (IOException e) {
          Utils.FileReadError(inputbin.getName());           }
   }
   
    private void ReadGridData(String file_lon, String file_lat){
        Pos p = new Pos();
        try{
        BufferedReader r = new BufferedReader(new FileReader(new File(file_lon)));
        String line = r.readLine();
        nlon = p.ReadInt(line);
        nlat = p.ReadInt(line);
        p.SetZero();
        for (int j=0; j<nlat; j++){
            line = r.readLine();
            for (int i=0; i<nlon; i++){
                xcrd[j][i] = p.ReadFloat(line);
            }
            p.SetZero();
        }
        }catch (IOException e) {
             System.out.println("Failed task...");    
             System.exit(0);
        }
        try{
        BufferedReader r = new BufferedReader(new FileReader(new File(file_lat)));
        String line = r.readLine();
        nlon = p.ReadInt(line);
        nlat = p.ReadInt(line);
        p.SetZero();        
        for (int j=nlat-1; j>=0; j--){
            line = r.readLine();
            for (int i=0; i<nlon; i++){
                ycrd[j][i] = p.ReadFloat(line);
            }
            p.SetZero();
        }
        }catch (IOException e) {
             System.out.println("Failed task...");    
             System.exit(0);
        }        
    }
    
    
    private void ReadDymFile(File ifile, int nbytetoskip){
       float[][] ivar = new float[nlon+2][nlat+2];
        ByteBuffer dim_buffer = null;
       //float mmin = 0;
       //float mmax = 0;
//       dim_buffer = ByteBuffer.allocate((nlon+2)*(nlat+2)*4);
       dim_buffer = ByteBuffer.allocate(nlon*nlat*4);
       try{
            FileInputStream fileInput = new FileInputStream(ifile);
            FileChannel fc = fileInput.getChannel();
            fc.position(nbytetoskip);
            fc.read(dim_buffer);
            dim_buffer.order(ByteOrder.LITTLE_ENDIAN);
            dim_buffer.rewind();
            
            //=== DATA =======================
            for (int j=0;j<nlat;j++){
                for (int i=0;i<nlon;i++){
                    var[i][j]= dim_buffer.getFloat();
                }
            }

/*            for (int j=0;j<nlat+2;j++){
                for (int i=0;i<nlon+2;i++){
                    ivar[i][j]= dim_buffer.getFloat();
                }
            }
            for (int j=0;j<nlat;j++){
                for (int i=0;i<nlon;i++){
                    var[i][j] = ivar[i+1][j+2];
                }
            }
  */          
            //mmin = Utils.getMIN(var, nlon, nlat);
            //mmax = Utils.getMAX(var, nlon, nlat);
        } catch (IOException e) {
            Utils.FileReadError(inputbin.getName());
           }
    }
    
    private void Store3d_var(float[][] data, int n){
      for (int i=0; i<size_x; i++)
           for(int j=0; j<size_y; j++){
               all_series[i][j][n] = data[i][j];
           }
    }
    
     private void WriteDymFile(double[][] data, boolean FileMode) {
        try{
            if (FileMode){
                FileOutputStream ofile = new FileOutputStream(outputbin);
                if (Old_Format) 
                    WriteOldDymFileHeader(ofile);
                else
                    WriteDymFileHeader(ofile);
                if (DoublePrecision)
                    WriteDataInDoubles(ofile, data);
                else
                    WriteDataInFloats(ofile, data);
                ofile.close();
            }
            else {
                FileOutputStream ofile = new FileOutputStream(outputbin,true);
                if (DoublePrecision)
                    WriteDataInDoubles(ofile, data);
                else
                    WriteDataInFloats(ofile, data);
                ofile.close();
            }
        } catch (IOException e) {
            Utils.FileWriteError(outputbin.getName());
            } 
    }
    
        private void WriteDymFile(float[][] data, boolean FileMode) {
        double[][] Data = new double[nlon][nlat];
        if (DoublePrecision){
        for (int i = 0; i < nlon; i++) 
            for (int j = 0; j < nlat; j++) {
                Data[i][j] = data[i][j];
            }
        }
        try{
            if (FileMode){
                FileOutputStream ofile = new FileOutputStream(outputbin);
                if (Old_Format) 
                    WriteOldDymFileHeader(ofile);
                else
                    WriteDymFileHeader(ofile);
                if (DoublePrecision)
                    WriteDataInDoubles(ofile, Data);
                else
                    WriteDataInFloats(ofile, data);
                ofile.close();
            }
            else {
                FileOutputStream ofile = new FileOutputStream(outputbin,true);
                if (DoublePrecision)
                    WriteDataInDoubles(ofile, Data);
                else
                    WriteDataInFloats(ofile, data);
                ofile.close();
            }
        } catch (IOException e) {
            Utils.FileWriteError(outputbin.getName());
            } 
    }
        
    private void WriteDymFileHeader(OutputStream out) throws IOException {
        
	//=== Magic ID, equal to 'DYM2' =========
        char[] ch = new char[4];
        for (int i=0; i<ch.length; i++)
            ch[i] = Magic_ID.charAt(i);
        Utils.writeChar(out,ch);       
       
	//=== Stretching function ID ============ 
	Utils.writeInt(out,FUNC);
	
        //=== MIN_VALUE =========================
        //minValue = 0.f;
        //minValue = Utils.getMIN(all_series,size_x,size_y,nbt_total);
        Utils.writeFloat(out,(float)minValue);
        
        //=== MAX_VALUE =========================
        //maxValue = 0.f;
        //maxValue = Utils.getMAX(all_series,size_x,size_y,nbt_total);
        Utils.writeFloat(out,(float)maxValue);
        
        //=== NLON ==============================
        Utils.writeInt(out,size_x);

        //=== NLAT ==============================
        Utils.writeInt(out,size_y);
        
        //=== NLEVEL ============================
        Utils.writeInt(out,nlevel);

        //=== START_DATE ========================
        Utils.writeFloat(out,startdate);
       
        //=== END_DATE ==========================
        Utils.writeFloat(out,enddate);

        //=== XLON ==============================
        float[] bufferF_x = new float[size_x];
        for (int i = 0; i < size_y; i++) {
            for (int j = 0; j < size_x; j++) {
                 bufferF_x[j] = (float)coord[j][i].x;//xcrd[i][j];
            }
        Utils.writeFloat(out,bufferF_x);
        }

        //=== YLAT ==============================
        float[] bufferF_y = new float[size_x];
        for (int i = 0; i < size_y; i++) {
            for (int j = 0; j < size_x; j++) {
                bufferF_y[j] = (float)coord[j][i].y;//ycrd[i][j];
            }
        Utils.writeFloat(out,bufferF_y);
        }
        //=== ZLEVEL ============================
        Utils.writeFloat(out,zlevel);

        //=== MASK ==============================
        int[] bufferI = new int[size_x];
        for (int j = 0; j < size_y; j++) {
            for (int i = 0; i < size_x; i++) {
                bufferI[i] = mask[i][j];//msk[i][j];//mskt[i][j];
            }
            Utils.writeInt(out,bufferI);
        }
    }
    
       private void WriteOldDymFileHeader(OutputStream out) throws IOException {
        //=== NLON ==============================
        Utils.writeInt(out,size_x);

        //=== NLAT ==============================
        Utils.writeInt(out,size_y);
        
        //=== NLEVEL ============================
        Utils.writeInt(out,nlevel);

        //=== START_DATE ========================
        Utils.writeFloat(out,startdate);
       
        //=== END_DATE ==========================
        Utils.writeFloat(out,enddate);

        //=== XLON ==============================
        float[] bufferF_x = new float[size_x];
        for (int i = 0; i < size_x; i++) {
             bufferF_x[i] = (float)coord[i][0].x;
        }
        Utils.writeFloat(out,bufferF_x);

        //=== YLAT ==============================
        float[] bufferF_y = new float[size_y];
        for (int j = 0; j < size_y; j++) {
             bufferF_y[j] = (float)coord[0][j].y;
        }
        Utils.writeFloat(out,bufferF_y);
        
        //=== ZLEVEL ============================
        Utils.writeFloat(out,zlevel);

        //=== MASK ==============================
        int[] bufferI = new int[size_x];
        for (int j = 0; j < size_y; j++) {
            for (int i = 0; i < size_x; i++) {
                bufferI[i] = mask[i][j];
            }
            Utils.writeInt(out,bufferI);
        }
    }
    
    private void WriteDataInFloats(FileOutputStream out, double[][] mat2d) throws IOException {
        //=== DATA ==========================
        float[] bufferF = new float[size_x];
            for (int j = 0; j < size_y; j++) {
                for (int i = 0; i < size_x; i++) {
                    bufferF[i] = (float)mat2d[i][j];
                }
                Utils.writeFloat(out,bufferF);
            }
    }
    
    private void WriteDataInFloats(FileOutputStream out, float[][] mat2d) throws IOException {
        //=== DATA ==========================
        float[] bufferF = new float[size_x];
            for (int j = 0; j < size_y; j++) {
                for (int i = 0; i < size_x; i++) {
                    bufferF[i] = (float)mat2d[i][j];
                }
                Utils.writeFloat(out,bufferF);
            }
    }

    private void WriteDataInDoubles(FileOutputStream out, double[][] mat2d) throws IOException {
        //=== DATA ==========================
        double[] bufferF = new double[size_x];
            for (int j = 0; j < size_y; j++) {
                for (int i = 0; i < size_x; i++) {
                    bufferF[i] = mat2d[i][j];
                }
                Utils.writeDouble(out,bufferF);
            }
    }
    
     private void ScaleData(float[][] var, int nx, int ny, float max){        
        for (int i=0; i<nx; i++){
            for (int j=0; j<ny; j++){
                var[i][j] = var[i][j]/max;
            }
        }
    }
    
    private void ScaleData(double[][] var, int nx, int ny, float min_init, float max_init){
        //double min_cur = Utils.getMIN(var, nx, ny, mask);
        for (int i=0; i<nx; i++){
            for (int j=0; j<ny; j++){
                if (mask[i][j]!=0){
                    if (var[i][j]>max_init)
                        var[i][j] = max_init;
                    if (var[i][j]<min_init)
                        var[i][j] = min_init;
                }
            }
        }
    }
    
    private void InitInterpolation(){
        int nmax;
	if (nlat < nlon) 
            nmax = nlon;
        else 
            nmax = nlat;
        xlon_im  = new float[nlon];
        ylat_im  = new float[nlat];
        for (int i=0; i<nlon; i++)
            xlon_im[i] = xlon[i]-xlon[0];
        for (int j=0; j<nlat; j++)
            ylat_im[j] = ylat[0]-ylat[j];
        vmid  = new double[nlon][nlat];
	vfin  = new double[size_x][size_y];
        y2    = new double[nmax];
	ytmp  = new double[nmax];
	y2tmp = new double[nmax];
	yytmp = new double[nmax];
	u     = new double[nmax];
        if (FUNC != 2 || Bilinear){
            Bisection();
        }
        if (AdjustMsk)
            AdjustMask();        
}
    private void AdjustMask(){
        int i,j,il,jl;
        for (i=0; i<size_x; i++){
            for (j=0; j<size_y; j++){
                if (mask[i][j] != 0){
                    il = xind[i];//index of the left coordinate on the initial grid 
                    jl = yind[j]; 
                    if (mskt[il][jl]==0 || mskt[il+1][jl]==0
                        || mskt[il+1][jl+1]==0 || mskt[il][jl+1]==0)
                        mask[i][j] = 0;
                }
            }                    
	}
    }
    
    private void Bisection(){
        //works for rectangular grids only
        xind  = new int[size_x];
        yind  = new int[size_y];
        float[] y;
        y = new float[nlat];
        int klo,khi,k;
        int i,j;
        for (i=0; i<size_x; i++){
            xind[i]=locate(xlon,nlon,coord[i][0].x);
        }
        for (j=0; j<size_y; j++){
            //since we have negative coordinates on south
            double crd_y = (coord[0][0].y-coord[0][j].y);         
            for (int jj=0; jj<nlat; jj++)
                y[jj] = (float)(ylat[0]-ylat[jj]);
            yind[j]=locate(y,nlat,crd_y);
        }
    }
    
    private int locate(float[] xx, int n, double x){
        int klo,khi,k;
        klo = 0;
        khi = n-1;
        while (khi-klo > 1) {
            k = (khi+klo) >> 1; //compute a midpoint
            if (x>xx[k]) 
                klo=k; 
            else 
                khi=k;
        }
        return klo;
    }

    private void Interpolate(float[][] var){
        if (Bilinear) {
            BilinearInt(var);
            return;
        } 
        if (Averaging){
            DoAveraging(var);
            return;
        }
        //=== let's calculate vmid
        Splie2(ylat_im,var,nlon,nlat);	
	if (FUNC == 2){ // if stretching is local we'll use metric...
            for (int i=0; i<size_x; i++){
		for (int j=0; j<size_y; j++){
                    if (mask[i][j] != 0){
                        if ((metric[i][j].x == 1)&&(metric[i][j].y == 1)){               
                            vfin[i][j] = var[i][j];
                        } else {
                            double crd_x = (coord[i][j].x-xlon[0]);
                            double crd_y = (ylat[0]-coord[i][j].y);
                            vfin[i][j] = Splin2(xlon_im,ylat_im,var,crd_x,crd_y,nlon,nlat);
                        }
                    }
                }
            }
        } else {
            for (int i=0; i<size_x; i++){
		for (int j=0; j<size_y; j++){
                    if (mask[i][j] != 0){
                     //=== get final interpolated data
                    //scale grid coordinates from 0 to nlon and nlat
                    //double crd_x = (coord[i][j].x-xlon[0]);
                    //double crd_y = (ylat[0]-coord[i][j].y);
                    double crd_x = (coord[i][j].x-coord[0][j].x);
                    double crd_y = (coord[i][0].y-coord[i][j].y);
                    //int jj = (int)Math.round(crd_y); 
                    Point p = get_indexes(coord[i][j].x,coord[i][j].y);
                    int jj=p.y;
                    //if (jj>148) jj = 148;
                    for (int ii=0; ii<nlon; ii++){
                        xlon_im[ii] = (float)(xcrd[jj][ii]-coord[0][j].x);
                        //if (xlon_im[ii]<0) xlon_im[ii] = 0;
                    }
                    //int ii = (int)Math.round(crd_x/2);
                    int ii=p.x; 
                    for (jj=0; jj<nlat; jj++){
                        ylat_im[jj] = (float)(coord[i][0].y-ycrd[jj][ii]);
                        //if (ylat_im[jj]<0) ylat_im[jj] = 0;
                    }                    
                    //vfin[i][j] = Splin2(xlon_im,ylat_im,var,crd_x,crd_y,nlon,nlat);
                    vfin[i][j] = Splin2_2(xlon_im,ylat_im,var,crd_x,crd_y,nlon,nlat,i,j);
                    }
                }
            }
        }
    }
    
    private Point get_indexes(double x, double y){
        Point indexes = new Point(0,0);
        for (int i=1; i<nlat; i++){
            for (int j=1; j<nlon; j++){        
                if ((xcrd[i][j-1]<=x)&&(x<=xcrd[i][j])
                    &&(ycrd[i-1][j]>=y)&&(y>=ycrd[i][j])){
                    indexes.x = j;
                    indexes.y = i;
                    break;
                }
            }
        }
        return indexes;
    }

    private double Splint(float[] crd, double[] ya, double[] y2a, int n, double x){
        //bisection: find the locations of unknown data points within the initial grid
        int klo,khi,k;
        double h,b,a;
        klo = 0;
        khi = n-1;
        while (khi-klo > 1) {
            k = Math.round((khi+klo)/2);
            if (crd[k] > x) 
                khi=k; 
            else 
                klo=k;
	}
        h=crd[khi]-crd[klo];
        a= (crd[khi]-x)/h;
        b= (x-crd[klo])/h;
        //now evaluate and return cubic spline polynomial
        return a*ya[klo]+b*ya[khi]+(a*(a*a-1)*y2a[klo]+b*(b*b-1)*y2a[khi])*h*h/6;
}
       private double Splint2(float[] crd, double[] ya, double[] y2a, int n, double x, int k){
        //bisection done before
        double h,a,b;
        h=crd[k+1]-crd[k];
        a= (crd[k+1]-x)/h;
        b= (x-crd[k])/h;
        //now evaluate and return cubic spline polynomial
        return a*ya[k]+b*ya[k+1]+(a*(a*a-1)*y2a[k]+b*(b*b-1)*y2a[k+1])*h*h/6;
}

    private double[] Spline(float[] x, double[] y, int n){
	int i,j;
	double p,qn,sig,un;
        
	y2[0]=0; 
	u[0] =0;  
	for(i=1; i<n-1; i++)
	{
        sig   = (x[i]-x[i-1])/(x[i+1]-x[i-1]);
        p     = sig*y2[i-1]+2;
        y2[i] = (sig-1)/p;
        u[i]  = (y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]);
        u[i]  = (6*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
	}
	qn = 0; 
	un = 0;
	y2[n-1] = (un-qn*u[n-2])/(qn*y2[n-2]+1);

	for (j=n-2; j>=0; j--)
        y2[j] = y2[j]*y2[j+1]+u[j];

	return y2;
        
    }

    private void Splie2(float[] yycrd, float[][]v_pr, int m, int n){
        int i,j;
	for (i=0; i<m; i++){
		for (j=0; j<n; j++){
			ytmp[j] = v_pr[i][j];
                        yycrd[j]= ycrd[0][i]-ycrd[j][i];
                }
		Spline(yycrd,ytmp,n);	// returns y2
		for (j=0; j<n; j++)
			vmid[i][j] = y2[j];
	}
}

    private double Splin2(float[] xcrd, float[] ycrd, float[][] v_pr, double xnew, double ynew, int m, int n){ //x1 and x2 are the new coordinates
        int i,j;
        for (i=0; i<m; i++){
		for (j=0; j<n; j++){
			ytmp[j] = v_pr[i][j]; 
			y2tmp[j]= vmid[i][j];
		}
		yytmp[i] = Splint(ycrd,ytmp,y2tmp,n,ynew);
	}
	Spline(xcrd,yytmp,m); // returns y2
	return Splint(xcrd,yytmp,y2,m,xnew);
    }
    
    private double Splin2_2(float[] xcrd, float[] ycrd, float[][] v_pr, double xnew, double ynew, int m, int n, int ii, int jj){ //x1 and x2 are the new coordinates
        int i,j;
        for (i=0; i<m; i++){
		for (j=0; j<n; j++){
			ytmp[j] = v_pr[i][j]; 
			y2tmp[j]= vmid[i][j];
		}
		yytmp[i] = Splint2(ycrd,ytmp,y2tmp,n,ynew,yind[jj]);
	}
	Spline(xcrd,yytmp,m); // returns y2
	return Splint2(xcrd,yytmp,y2,m,xnew,xind[ii]);
    }
    
    private void BilinearInt(float[][] var){
        int i,j,il,jl;
        double x,y,f1,f2,f3,f4;
        for (i=0; i<size_x; i++){
            for (j=0; j<size_y; j++){
                if (mask[i][j] != 0){
                    il = xind[i];//index of the left coordinate on the initial grid 
                    jl = yind[j];                     
                    x = (coord[i][j].x-xlon[il])/(xlon[il+1]-xlon[il]);
                    y = (ylat[jl]-coord[i][j].y)/(ylat[jl]-ylat[jl+1]);
                    f1 = var[il][jl];     f2 = var[il+1][jl];
                    f3 = var[il+1][jl+1]; f4 = var[il][jl+1];
                    if (f1!=misval && f2!=misval && f3!=misval && f4!=misval)
                        vfin[i][j] = (1-x)*(1-y)*f1+x*(1-y)*f2+x*y*f3+(1-x)*y*f4;
                }
            }                    
	}
    }
    
    private void DoAveraging(float[][] var){
        int i,j,ii,jj,count;
        double x,y,sum;
        int m = nlon/size_x;
        int n = nlat/size_y;
        for (i=0; i<size_x; i++){
            for (j=0; j<size_y; j++){
                if (mask[i][j] != 0){
                    sum = 0;
                    count = 0;
                    for (ii=i*m; ii<(i+1)*m; ii++)
                        for (jj=j*n; jj<(j+1)*n; jj++)
                            if ((msk[jj][ii] != 0) || (var[ii][jj] != misval)){
                                sum += var[ii][jj]; 
                                count++;
                            }
                    if (count>0)
                        vfin[i][j] = sum/count;
                    else{
                        switch (MVmethod) {
                            case 0: 
                                vfin[i][j] = misval;   
                                break;
                            case 1:
                                vfin[i][j] = 0;
                                break;
                            case 2: 
                                GetInterpolatedValue();
                        }
                    }
                }
            }                    
	}        
    }
    
    private void GetInterpolatedValue(){
        //add code
    }
                
}

   
